#include "MatrixMultiply.h"
#include "Matrix.h"
#include <climits>
#include <ctime>
#include <cfloat>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <omp.h>
#include <math.h>
using namespace std;

extern unsigned int matrices;

int main(void)
{
	int maxCount = 800;
	double timeTableNaive[maxCount];
	double timeTableDDOT[maxCount];
	double timeTableDAXPY[maxCount];
	double timeTableDGEMV[maxCount];
	double timeTableDGEMM[maxCount];
//Multiplikation Naiv
	for(int n=1; n < maxCount; n++)	
	{
		int ld = (n%2) ? n : n+1;
		Matrix matA(n, ld);
		Matrix matB(n, ld);
		Matrix matC(n, ld);

		MatrixMultiply mm(&matA, &matB, &matC);
// 9.4412e-08*n*n*n + (-6.2391e-05*n*n) + 1.2162e-02*n -4.8020e-01
		int iter = max(1, (int)round(9.4412e-08*n*n*n + (-6.2391e-05*n*n) + 1.2162e-02*n -4.8020e-01));
		double time = (double)FLT_MAX;
		for(int rep = 0; rep < 4; rep++)
		{
			double start = omp_get_wtime();
			for(int its=0; its < iter; its++)			
				mm.calcNaive();

			double stop = omp_get_wtime();
			double t = ((double)(stop-start))/iter;
			time = min(time, t);
		}
		cout << n << " / " << matrices << endl;		
		timeTableNaive[n-1] = time;
	}
//Multiplikation DDOT
	for(int n=1; n < maxCount; n++)	
	{
		int ld = (n%2) ? n : n+1;
		Matrix matA(n, ld);
		Matrix matB(n, ld);
		Matrix matC(n, ld);

		MatrixMultiply mm(&matA, &matB, &matC);
//8.4482e-08*n*n*n -6.0386e-05*n*n + 1.1654e-02*n  -4.5749e-01
		int iter = max(1, (int)round(8.4482e-08*n*n*n -6.0386e-05*n*n + 1.1654e-02*n  -4.5749e-01));
		double time = (double)FLT_MAX;
		for(int rep = 0; rep < 4; rep++)
		{
			double start = omp_get_wtime();
			for(int its=0; its < iter; its++)			
				mm.cblas_myDdot();

			double stop = omp_get_wtime();
			double t = ((double)(stop-start))/iter;
			time = min(time, t);
		}
		cout << n << " / " << matrices << endl;		
		timeTableDDOT[n-1] = time;
	}
//Multiplikation DAXPY
	for(int n=1; n < maxCount; n++)	
	{
		int ld = (n%2) ? n : n+1;
		Matrix matA(n, ld);
		Matrix matB(n, ld);
		Matrix matC(n, ld);

		MatrixMultiply mm(&matA, &matB, &matC);
//5.5757e-09*n*n*n  -3.5435e-06*n*n + 7.9843e-04*n  -3.6451e-02
		int iter = max(1, (int)round(5.5757e-09*n*n*n  -3.5435e-06*n*n + 7.9843e-04*n  -3.6451e-02));
		double time = (double)FLT_MAX;
		for(int rep = 0; rep < 4; rep++)
		{
			double start = omp_get_wtime();
			for(int its=0; its < iter; its++)			
				mm.cblas_myDaxpy();

			double stop = omp_get_wtime();
			double t = ((double)(stop-start))/iter;
			time = min(time, t);
		}
		cout << n << " / " << matrices << endl;		
		timeTableDAXPY[n-1] = time;
	}

//Multiplikation DGEMV
	for(int n=1; n < maxCount; n++)	
	{
		int ld = (n%2) ? n : n+1;
		Matrix matA(n, ld);
		Matrix matB(n, ld);
		Matrix matC(n, ld);

		MatrixMultiply mm(&matA, &matB, &matC);
//  6.3698e-09*n*n*n  -4.1777e-06*n*n +  9.3187e-04*n  -4.2544e-02

		int iter = max(1, (int)round(6.3698e-09*n*n*n  -4.1777e-06*n*n +  9.3187e-04*n  -4.2544e-02));
		double time = (double)FLT_MAX;
		for(int rep = 0; rep < 4; rep++)
		{
			double start = omp_get_wtime();
			for(int its=0; its < iter; its++)			
				mm.cblas_myDGEMV();

			double stop = omp_get_wtime();
			double t = ((double)(stop-start))/iter;
			time = min(time, t);
		}
		cout << n << " / " << matrices << endl;		
		timeTableDGEMV[n-1] = time;
	}
//Multiplikation DGEMM
	for(int n=1; n < maxCount; n++)	
	{
		int ld = (n%2) ? n : n+1;
		Matrix matA(n, ld);
		Matrix matB(n, ld);
		Matrix matC(n, ld);

		MatrixMultiply mm(&matA, &matB, &matC);
// 2.4218e-10*n*n*n + 2.2926e-08*n*n -1.4242e-06*n + 1.2534e-04

		int iter = max(1, (int)round(2.4218e-10*n*n*n + 2.2926e-08*n*n -1.4242e-06*n + 1.2534e-04));
		double time = (double)FLT_MAX;
		for(int rep = 0; rep < 4; rep++)
		{
			double start = omp_get_wtime();
			for(int its=0; its < iter; its++)			
				mm.cblas_myDGEMM();

			double stop = omp_get_wtime();
			double t = ((double)(stop-start))/iter;
			time = min(time, t);
		}
		cout << n << " / " << matrices << endl;		
		timeTableDGEMM[n-1] = time;
	}
//	double *matCContent = matC->GetMatrix();
//	int sizeOfC = matC->GetLD()*matC->GetRowCount();
	fstream myfile;
	myfile.open("Rene Eichert ZeitenMessung.txt", ios::out | ios::binary | ios::trunc);
	//FileHead
	myfile << "Problemgroesse\tNaiv\tDDOT\tDAXPY\tDGEMV\tDGEMM" << endl;
	for(int i = 0; i < maxCount; i++)
	{
		myfile << i << "\t" << timeTableNaive[i] << "\t" << timeTableDDOT[i] << "\t" << timeTableDAXPY[i] << "\t" << timeTableDGEMV[i] << "\t" << timeTableDGEMM[i] << endl;
	}			
	myfile.close();
	cout << "Done with the Measuring, press Enter to Exit" << endl;
	cin.get();
	return 0;
}
